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C.  I 


ABSTRACT 

The  theory  and  background  of  the  algebraic  substitution 
synthesis  method  for  the  digital  filters  from  continuous 
filter  characterizations  is  presented  with  emphasis  on  the 
frequency  distortion  phenomenon.   An  analysis  of  the  Forward 
Euler,  Backward  Euler  and  Trapezoidal  numerical  integration 
algorithms  is  undertaken  and  appropriate  transformations 
are  obtained.   A  general  integration  formula,  encompassing 
the  above  algorithms  as  special  cases,  is  analyzed  and  its 
application  to  the  synthesis  problem  is  pointed  out.   Direct 
transformations  for  discrete  filter  frequency  response 
characteristics  from  continuous  filter  characterizations  are 
derived. 
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I.   INTRODUCTION 

Digital  filter  design  can  be  accomplished  using  several 
different  techniques.  [1]   One  particularly  convenient  method 
is  to  transform  a  continuous  filter  into  a  digital  filter  by 
means  of  some  substitution  procedure.   This  offers  the  designer 
a  large  body  of  filter  designs  from  which  to  draw  the  form 
applicable  to  a  particular  task.   However,  this  procedure  has 
certain  drawbacks.   Chief  among  these  is  the  distortion  which 
results  when  transforming  frequencies  from  the  continuous 
domain  into  frequencies  in  the  discrete  domain.   An  under- 
standing of  this  frequency  distortion  is  essential  when  the 
substitution  method  is  employed.   Fig.  1.1  shows  the  relation- 
ship between  the  continuous  filter  and  the  digital  filter. 

The  transfer  function  of  the  continuous  filter  is  derived 
from  a  differential  equation  by  means  of  the  Laplace  Transform. 
The  frequency  response  of  the  continuous  filter  can  be  eval- 
uated by  letting  s  =  joo  as  shown  in  the  upper  right  of  the 
diagram.   A  difference  equation  can  be  formed  by  applying  a 
numerical  integration  algorithm  to  the  continuous  differential 
equation.   This  difference  equation  can  then  be  z-transformed 
to  arrive  at  a  discrete  transfer  function  as  shown.   This 
discrete  transfer  function  can  also  be  developed  directly  by 
applying  a  substitution  (s=  fCz))  to  the  continuous  transfer 
function.   This  substitution  is  derived  from  the  numerical 
integration  algorithm  and  the  two  discrete  transfer  functions 
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thus  obtained  are  identical.   This  procedure  is  the  key  to  the 
algebraic  substitution  technique. 

The  frequency  response  of  the  discrete  transfer  function 
can  now  be  evaluated  by  letting  z  =  exp  ( j  £T)  .   The  spectrum 
thus  obtained  will  not,  in  general,  be  identical  to  that  of 
the  continuous  filter  spectrum.   The  difference  between  the  two 
gives  rise  to  frequency  distortion.   The  purpose  of  this  thesis 
is  to  analyze  this  frequency  distortion  for  several  numerical 
integration  algorithms . 

The  source  of  the  frequency  distortion  can  be  seen  in 
Fig.  1.2.   This  is  a  representation  of  the  effect  of  the  s-plane 
to  z-plane  transformation  on  the  sinusoidal  steady-state  fre- 
quency domain  characteristics  of  the  filter.   If  the  transfor- 
mation was  linear  as  depicted  by  the  dotted  line,  the  continuous 
spectrum  would  be  perfectly  reproduced  in  the  discrete  fre-' 
quency  domain.   However,  the  transformation  is  not  linear.   Its 
shape,  of  course,  depends  on  the  particular  integration  algorithm 
used  and  the  transforming  function  is  periodic  in  nature.   The 
solid  line  represents  an  arbitrary  transforming  function  and 
its  effect  on  the  frequency  characteristics  of  the  discrete 
filter.   The  distortion  that  results  from  this  nonlinear  trans- 
formation is  the  subject  of  the  analysis  that  follows. 

Much  use  is  made  of  the  Z  -  Transform  technique  to  accomplish 
this  analysis  and  the  theory  and  use  of  this  technique  is 
described  by  many  authors  including  E.  I.  Jury  [2]  and  B.  C. 
Kuo  [3].   A  paper  comparing  various  digital  integrators  was 
written  by  J.  R.  Salzer  [4]  and  is  one  of  the  classic  references 
on  this  subject.   Gold  and  Rader  [5]  have  discussed  the 
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Fig.  1.2.   GENERAL  TRANSFORMATION  OF  CONTINUOUS 

FREQUENCY  SPECTRUM  TO  DISCRETE  FREQUENCY 
SPECTRUM. 


substitution  technique  in  a  limited  form,  and  Golden  and 
Kaiser  [6]  have  discussed  the  bilinear  substitution,  which 
corresponds  to  trapezoidal  integration,  in  a  thorough  manner. 

Chapter  II  of  this  thesis  contains  the  theory  and  develop- 
ment of  the  algebraic  substitution  technique  and  describes  how 
a  filter  design  can  be  accomplished  using  this  method.   Chapter 
III  presents  a  detailed  analysis  of  three  specific  numerical 
integration  algorithms  -  Backward  and  Forward  Euler  and  the 
Trapezoidal  Rule.   The  effect  of  frequency  distortion  is  shown 
and  analyzed.   Also  specific  limitations  on  the  use  of  the 
Euler  methods  are  pointed  out.   Chapter  IV  discusses  a  general 
integration  formula  (which  includes  the  above  as  special  cases) 
and  how  it  can  be  employed.   General  case  sinusoidal  steady- 
state  frequency  transformations  are  then  derived,  and  an 
example  of  a  filter  synthesis  is  presented  to  demonstrate  the 
distortion  problem.   Chapter  V  contains  the  summary  and 
suggestions  for  further  research. 
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II.   THEORY  AND  DEVELOPMENT 
OF  ALGEBRAIC  SUBSTITUTION  DESIGN  METHOD 

There  are  many  instances  when  a  discrete  time  representa- 
tion of  a  continuous  system  or  filter  is  required.   Digital 
simulation  of  a  control  system  for  computer  analysis  is  one 
example.   Another  important  example  is  the  extension  of  con- 
tinuous filter  theory  to  the  design  and  synthesis  of  digital 
filters.   This  application  in  particular   provides  incentive 
for  developing  a  method  of  design  which  will  incorporate  the 
great  wealth  of  information  that  has  been  accumulated  for 
continuous  filters.   Designs  for  continuous  filters  are  well- 
known  and  have  been  tabulated  in  various  handbooks  [7] .   Thus 
a  direct  method  for  obtaining  a  digital  realization  of  a 
continuous  filter  is  not  only  theoretically  appealing  but  has 
definite  practical  benefits. 

The  usual  procedure  for  implementing  a  digital  filter  is 
to  first  express  the  filter  as  a  transfer  function  in  the  z- 
domain.   This  transfer  function  is  normally  expressed  as  a 
ratio  of  two  polynomials  in  inverse  powers  of  z,  where  z 
represents  a  time  delay  [8].   This  ratio  then  indicates  the 
linear  operations  which  must  be  performed  upon  the  past  values 
of  the  input  and  output  samples  to  obtain  the  present  value 
of  the  output.   For  example,  consider  the ' following  transfer 
function 

H(Z)  =  |gf  (2-1) 
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Since  this  transfer  function  is  the  ratio  of  the  transform 
of  the  output  to  the  transform  of  the  input  this  can  be 
expressed  as 

Z      a.  z 

Y  r~\         A(z)      i=Q   X ,0  -. 

X  (z)  -  BTiT n 7  (2"2) 

1  +  Z      b.  z"D 

j-l   ^ 

The  latter  expression  is  the  recursive  form  of  the  digital 
filter  transfer  function. 

In  a  like  manner,  the  most  common  description  of  continuous 
filters  in  the  Laplace  or  s-domain  is  as  the  ratio  of  two 
polynomials  in  s.   Consequently  any  substitution  that  can  be 
used  to  directly  replace  the  variable  s  in  the  continuous 
transfer  function  with  a  function  of  z  enables  a  discrete 
filter  realization  to  be  formed  directly.   Thus  if 

s  =  f(z)  (2-3) 

Then  H(s)  becomes 


G(z)  =  H(s)  =  H[f(z)]  (2-4) 

s  =  f(z) 


Where  H(s)  is  the  continuous  filter  transfer  function  and 
f(z)  is  a  generic  algebraic  expression. 

If  the  function  f Cz)  is  a  ratio  of  two  polynomials  of 
inverse  powers  of  z,  and  if  H(s)  is  a  ratio  of  two  polynomials 
in  s,  then  it  is  a  simple  matter  to  show  that  the  G(z)  thus 
obtained  will  always  be  the  ratio  of  two  polynomials  of 
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inverse  powers  of  z.   The  process  expressed  by  Eq.  (2-4)  is  an 
algebraic  substitution  method  which  converts  a  rational 
polynomial  in  s  to  a  rational  polynomial  in  z.   It  is  the 
nature  of  this  direct  substitution  which  is  discussed  in  detail. 

A.   DEVELOPMENT  OF  THE  ALGEBRAIC  SUBSTITUTION 

The  question  that  must  now  be  answered  is  what  substitu- 
tions can  be  used  in  Eq.  (2-4)  and  how  can  they  be  obtained? 
To  answer  this  question  properly,  it  is  necessary  to  examine 
first  the  nature  of  the  Laplace  variable  s.   This  variable 
represents  an  operation  in  the  time  domain.   S  is  a  differentia- 
tion and   1/s  is  an  integration.   Thus  it  is  intuitively  felt 
that  whatever  f (z)  is  used,  it  should  perform  an  approximate 
operation  in  the  discrete  time  domain  similar  to  that  of  the 
Laplace  operator  in  the  continuous  time  domain.   Partarrieu 
[9]  has  shown  that  the  selection  of  a  particular  function  of 
Z  to  substitute  for  1/s  corresponds  to  the  adoption  of  a 
discrete  time  numerical  integration  algorithm.   A  brief  review 
of  his  development  will  be  helpful. 

Consider  a  continuous  transfer  function 

an  +  a,  s  +  .  .  .  +  a   s 

H(s)  =  — — -  ;n  >  m  (2-5) 

b„  +  bn  s  +  . . .  +  b   s 
0     1  n 

This  transfer  function  represents  the  ratio  of  the  transform 
of  the  output,  Y(s),  to  the  transform  of  the  input  signal,  X(s) 


H(s)  =  |£ff  (2-6) 
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This  equation  can  be  expressed  as  follows 


Y(s)  _  W(s)  ,  YCs)  _ 


where  W(s)  is  a  signal  defined  by  the  following  equations 


W(s)   .  1 

X(s)    b.  +  b.  s  +  ...  +  b   sn 
0     1  n 


(2-8) 


Y(s)  _  _   L  .   „  _,      ■  i   -m 

m 


Wtflt  =  ao  +  ai  s  +  *••  +  a~  s  (2"9) 


Rearranging  Eqs .  (2-8)  and  (2-9)  gives 


n-1 
b   sn  W(s)  =  X(s)  -  I      b.  s1  W(s)  (2-10) 

i=0 


m 
Y(s)  =  Z      a.  sD  W(s)  (2-11) 

j=0   => 


Equations  (2-10)  and  (2-11)  can  be  depicted  in  an  analog 
computer  simulation  as  shown  in  Fig.  2.1.   It  can  be  seen  that 
the  replacement  of  the  s    factors  (integrators)  by  a  function 
of  z  denoted  by  [f(z)]    corresponds  to  a  direct  algebraic 
substitution  for  the  variable  s.   If  the  [f(z)]    is  selected 
from  the  many  numerical  integration  algorithms  available  in 
the  literature  [10] ,  it  follows  that  the  substitution  will,  in 
fact,  correspond  to  the  operation  of  integration.   How  good 
a  job  it  does,  in  the  sense  as  to  how  well  the  discrete  time 
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approximation  represents  the  continuous  transfer  function, 
depends  upon  which,  integration  algorithm  was  selected.   A 
study  of  several  possible  integration  schemes  and  the  effect 
of  their  approximations  in  the  sinusoidal  steady-state  fre- 
quency domain  forms  the  basis  of  this  thesis. 

B.   FREQUENCY  RESPONSE  CONSIDERATIONS 

A  problem  that  the  filter  designer  must  face  when  formulat- 
ing a  filter  transfer  function  by  the  algebraic  substitution 
method  is  that  of  translation  from  the  S  and  the  Z  domain  to 
the  sinusoidal  steady-state  frequency  domains.   A  common  way 
to  obtain  the  frequency  response  of  a  continuous  filter  is  to 
let 

s  =  j  to  (2-12) 

and  evaluate  the  magnitude  squared  and  phase  characteristics 
of  the  transfer  function.   For  example,  suppose  the  desired 
transfer  function  is  as  shown 

H(s)  =  |{fl  (2-13) 

where  A(s)  and  B(s)  are  rational  polynomials  in  s,  then  the 
magnitude  function  would  be 

lH(jco)|2  =  \*W\l    =   M^  (2-14) 

|B(jw)  \Z         B(oj  ) 

This  expression  can  then  be  plotted  and  the  filter  classified 
as  one  of  several  types  (low  pass,  band  pass,  band  stop,  etc.) 
An  example  of  this  procedure  is  shown  in  Fig.  2.2. 
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Fig.  2.2.   EXAMPLE  OF  FORMULATING  FREQUENCY 
CHARACTERISTICS  FROM  CONTINUOUS 
TRANSFER  FUNCTION. 
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A  similar  concept  of  filter  frequency  response  exists  for 
the  digital  filter.   A  designer  will  want  to  be  able  to  examine 
the  response  from  the  specific  filter  and  determine  if,  in  the 
frequency  domain,  it  meets  the  needs  for  which  it  has  been 
designed.   To  do  this  it  is  necessary  to  recall  that  [2] 


z  =  exp  (sT)   so  that   z  =  ejQT  (2-15) 


where  T  is  the  sampling  period  and  s  =  jfi.   The  symbol  0,    is  used 
here  rather  than  co  so  as  to  distinguish  the  sinusoidal  steady- 
state  frequency  for  the  continuous  transfer  function  from  the 
sinusoidal  steady-state  frequency  for  the  discrete  time  trans- 
fer function.   Thus,  to  facilitate  the  formulation  of  the 
frequency  response  characteristics  of  the  discrete  transfer 
function  directly  from  the  sinusoidal  steady-state  character- 
istics of  the  continuous  transfer  function  it  is  convenient 
to  have  a  second  substitution 

to  =  F(jft)  (2-16) 

that  can  be  inserted  in  Eq.  (2-14)  to  arrive  directly  at  the 
magnitude  squared  frequency  domain  transfer  function  for  the 
discrete  case.  A  block  diagram  of  this  procedure  is  shown  in 
Fig.  2.3.  The  derivation  of  these  frequency-related  functions 
along  with  the  resulting  distortion  of  the  filter  character- 
istics of  the  transfer  function  is  discussed  in  detail  later 
on. 

In  the  process  of  translating  frequencies  from  one  domain 
to  another  the  designer  encounters  the  question  of  whether 
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these  frequency  translations  will  be  linear.   In  general  they 
will  not  be  linear  and  therefore  a  certain  amount  of  distortion 
or  "warping"  of  the  frequency  scales  will  result.   This 
phenomenon  poses  certain  restrictions  and  limitations  on  the 
use  of  the  algebraic  substitution  as  will  be  discussed. 

C.   DIRECT  PROCEDURE  FOR  FILTER  SYNTHESIS 

The  object  of  this  section  is  to  provide  a  detailed  outline 
of  the  steps  to  be  followed  in  employing  the  algebraic  substi- 
tution method  for  the  direct  synthesis  of  a  digital  filter. 

1.  Determine  the  specifications. 

Before  any  actual  design  work  can  be  done  the  designer 
must  thoroughly  understand  his  particular  problem  and  ascertain 
what  functions  his  filter  must  fulfill.   The  actual  specifica- 
tions should  be  carefully  studied  and  tabulated.   These  speci- 
fications may  include  cut-off  frequencies,  center  frequencies, 
roll-off  per  octave,  etc.   The  designer  must  have  a  clear 
picture  of  what  his  filter  must  do. 

2.  Solve  the  approximation  problem. 

This  step  encompasses  the  selection  of  the  most  suitable 
continuous  filter  that  will  approximate  the  desired  filter. 
There  are  a  wide  variety  of  filters  such  as  Butterworth, 
Chebyshev,  Elliptic,  and  others  that  are  described  in  detail 
and  their  characteristics  tabulated  in  various  filter  hand- 
books.  The  designer  must  be  able  to  choose  one  of  these 
classical  filter  designs  to  suit  his  particular  problem  and 
consequently  he  must  be  familiar  with  their  characteristics. 
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3.  Determine  the  transfer  function. 

The  transfer  function  of  the  selected  continuous  filter 
must  be  expressed  as  a  rational  function  of  s  in  order  to  apply 

the  algebraic  substitution  method.   In  addition  the  transfer 

2 

function  should  also  be  expressed  as  a  function  of  go  (or  go  ) 

in  order  to  have  an  insight  into  the  frequency  response  of  the 
filter. 

4.  Choose  the  numerical  integration  algorithm  for  the 
replacement  of  s. 

This  step  is  probably  the  most  important  one  in  the  chain. 
Whether  the  filter  performs  as  it  should  depends  to  a  great 
extent  on  the  proper  selection  of  the  numerical  integration 
algorithm.   A  familiarity  with  the  various  algorithms  available 
for  this  purpose  is  essential. 

5.  Perform  the  required  algebra. 

Once  the  required  substitution  functions  are  known  it  is  a 
relatively  easy  matter  to  accomplish  the  algebra  required  to 
put  the  new  G(z)  in  the  form  of  the  ratio  of  two  polynomials 
in  inverse  powers  of  z.   The  designer  may  have  to  make 
allowances  for  the  frequency  distortion  resulting  from  the 
non-linear  frequency  scale  translation.   The  required  frequency 
function  G(jfi)  can  also  be  formed.   This  function  will  in 
general  be  fairly  complex  and  will  require  effort  for  inter- 
pretation. 

6.  Implement  the  transfer  function. 

Once  the  necessary  algebra  has  been  accomplished  and  the 
transfer  function  has  been  formed  it  must  be  synthesized. 
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There  are  many  references  concerning  the  procedures  to  be 
employed  in  actually  synthesizing  the  filter  [1,  3] .   The 
actual  mechanics  of  this  process  are  not  within  the  purview 
of  this  thesis.   Hess  [11]  has  enumerated  the  twenty-four 
canonical  forms  for  the  synthesis  of  second-order  sections 
and  these  may  be  either  cascaded  or  paralleled  to  achieve  the 
necessary  implementation. 
7.   Test  the  completed  filter  design. 

One  method  of  observing  the  performance  of  the  designed 
filter  is  to  simulate  the  G(z)  on  a  digital  computer.   Programs 
are  available  for  this  purpose. 

It  can  be  seen  from  the  foregoing  that  the  critical  points 
in  the  design  procedure  are  incorporated  in  steps  4  and  5. 
The  remainder  of  this  thesis  is  directed  to  the  study  of 
integration  algorithms  and  their  application  to  the  design 
process . 
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III.   NUMERICAL  INTEGRATION  - 
SPECIFIC  ALGORITHMS 

Numerical  integration  is  the  process  of  computing  the  value 
of  a  definite  integral  from  a  set  of  numerical  values  of  the 
integrand.   This  process  is  accomplished  by  representing  the 
integrand  by  some  interpolation  formula  and  then  integrating 
this  formula  between  the  desired  limits.   The  integrand  is 
usually  represented  by  a  polynomial  which  approximates  the 
actual  shape  of  the  function  over  the  desired  interval.   Nearly 
all  of  the  numerical  integration  methods  employ  this  approx- 
imating technique  and  the  basic  difference  between  them  is  the 
complexity  of  the  approximations.   While,  in  general,  the  more 
complex  an  algorithm  is,  the  more  accurate  it  is;  this  does 
not  mean  that  it  is  more  desirable  to  use.   This  greater 
complexity  also  means  more  arithmetic  operations  to  be  performed 
and  an  increase  in  computation  time.   These  may  make  a  particu- 
lar algorithm  unsuitable. 

This  chapter  is  a  study  of  three  specific  algorithms  — 
Euler  methods,  forward  and  backward,  and  the  Trapezoidal 
method.   All  three  can  be  expressed  by  the  following  discrete 
integration  formula  [12] 

Y  "  J      fCt)  dt  (3-1) 

T 

k-1 


Yk  =  Yk-1  +  T  X    f  (Tk_]_)  +  T(l-A)  f  (Tk)  (3-2) 
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where  Y,  and  Y,  ,  are  the  computed  values  of  the  integral  at 
time  T,  and  time  T-  ,,  respectively;  fCT,)  is  the  value  of  the 
integrand  at  time  T,  .   A  is  the  parameter  that  determines  the 
particular  integration  scheme.   Consider  Eq.  (3-2)  with  the 
values  of  A  =  0,  1,  y,  respectively.   When 

A  =  0  ;  Yk  =  Yk_x  +  T  f (T,)     Backward  Euler        (3-3) 
A  =  1  ;  Yk  =  Yk-1  +  T  fCTk_1)   Forward  Euler         (3-4) 

A  =  j;Yk  =  Yk_x  +  |  If(Tk)  +  fCTj^)]  Trapezoid    (3-5) 

These  are  the  three  algorithms  and  it  can  be  seen  that  they 
differ  only  in  how  the  input  f  (t)  is  handled.   Fig.  3.1  is  a 
representation  of  how  these  three  algorithms  approximate  a 
function  for  integration.   Each  of  the  three  methods  will  be 
discussed  in  the  ensuing  sections. 

A.   FORWARD  EULER  INTEGRATION  METHOD 

The  Euler  methods  (both  forward  and  backward)  are  also 
known  as  the  rectangular  rule.   This  stems  from  the  geometric 
interpretation  of  the  integration  scheme.   The  integral  of,  or 
the  area  under,  a  curve  can  be  approximated  by   rectangles  as 
shown  in  Fig.  3.1.   The  area  is  brought  up-to-date  by  adding 
the   rectangle  indicated  to  the  area  already  accumulated. 
Thus  from  Eq.  (3-4) 

\  =  Yk-i  +  T  f (Tk-i>  (3'6) 


where  f  (T.  , )  is  the  value  of  the  function  being  integrated  at 

time  T.  ..  .   The  backward  Euler  as  discussed  in  the  next  section 
k-1 
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A        B  t 

0;Yk-Yk_1  =  Area  A  G  H  B  Backward  Euler 

2;Yk~Yk-l  =  Area  A  E  F  B  Trapezoidal 

1;Yk~Yk-l  =  Area  A  C  D  B  Forward  Euler 


Fig.  3.1.   GENERAL  INTEGRATION  FORMULA  FOR 
SPECIFIC  ALGORITHMS. 
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differs  from  the  forward  by  using  the  present  value  of  the 
input  f(Tj.)  to  compute  the  present  value  of  the  output  Y,  . 
It  can  be  seen  that  using  the  Euler  method  corresponds  to 
approximating  the  input  function  by  one  that  is  piecewise 
constant.   As  a  result  the  method  is  not  very  accurate  when 
the  function  is  changing  very  rapidly.   In  addition  the 
integration  time  step  will  have  to  be  small  for  the  approxima- 
tion  to  have  any  validity. 

To  derive  the  algebraic  substitution  for  the  forward  Euler 
integration  scheme,  it  is  necessary  to  represent  the  discrete 
difference  equation  (3-6)  in  the  z-domain.   Recalling  that 

z"1  =  exp  (-sT)  (3-7) 

and  that 

Z  [f(n-a)T]  =  z"a   Z[f(nT)]  (3-8) 


it  is  possible  to  z-transform  Eq .  (3-6)  yielding 

Y(z)  =  z"1  Y(z)  +  z"1  T  F(z)  (3-9) 


The  integration  transfer  function  is  then 

I    (z)  =  5L_z  (3-10) 

F       1-z"1 

This  transfer  function  representing  integration  in  the  s- 
domain  is 

\   (s)  -  i  '  (3-1D 
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A  comparison  of  Eqs .  (3-10)  and  (3-11)  indicates  the  following 
correspondence  can  be  made. 

1-z-1 
s  ■+  ±-Ar-  (3-12) 

Tz 

This  then  becomes  the  algebraic  substitution  for  s  for  the 
forward  Euler  integration  algorithm. 

It  is  now  possible  to  examine  the  sinusoidal  steady-state 
frequency  domain  characteristics  of  this  s-plane  to  z-plane 
transformation.   Letting  s  =  jfi  in  Eq.  (3-7)  and  applying  this 
to  Eq.  (3-12)  leads  to 

i  _  ^-J^T 
F(jfl)  =     _%T  (3-13) 

=  ^-  eJ      sin  2~  (3-14) 


F(jfi)  =  ^  [sin  QT  +   j (1-cos  OT)]  (3-15) 


This  then  is  the  required  F(jft)  for  the  forward  Euler.   This 
can  be  compared  with  the  ideal  integrator  in  the  continuous 
sense  to  discover  the  frequency  distortion  which  results.   The 
magnitude  and  phase  of  Eq.  (3-15)  are 

|F(jfi) |  =  |  sin  y~  (3-16) 

ARG  F(jft)  =  5^  +  |  (3-17) 

These   expressions  can  be  compared  to  the  magnitude  and  phase 
of  the  Laplace  variable  s  when  s  =  jco. 


27 


H(s)  =  s  (3-18) 

.   H(ju))  =  joo  (3-19) 

|HCju)  |  =  w  (3-20) 

ARG  H(joo)  =  j  (3-21) 

Eqs .  (3-16)  and  (3-20)  can  be  combined  to  provide  an  analysis 
of  the  magnitude  distortion  between  the  continuous  and  discrete 
operations . 


o)T     .   ftT  ,_  0_, 

2~  =  sin  y~  (3-22) 


This  deviation  from  linearity  for  the  magnitude  as  well  as  the 
phase  angle  deviation  is  shown  graphically  in  Fig.  3.2. 
Although  Eq.  (3-22)  is  nonlinear  it  can  be  seen  from  Fig.  3.2a 
that  there  is  a  region  of  linearity  close  to  the  origin.   Thus 
for  small  values  of  radian  frequency  the  forward  Euler  magni- 
tude spectrum  will  closely  approximate  that  of  the  ideal  in- 
tegrator.  This  can  be  made  more  quantitative  by  considering 
the  series  expansion  of  Eq.  (3-16) 

|F(jfl)|  =  %     [jr-  -  ^jf-  +  ...]  (3-23) 

3  3 

ftT 
a  o  -  u    x 

"     24 

Therefore  the  magnitude  spectrum  of  the  forward  Euler  will 

closely  resemble  that  of  the  ideal  integration  for  values  of 

ft  such  that 

3  2 
A  »  £_£_  (3-25) 

24 
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*F(jfi) 


ideal  integrator 


T 


a)  magnitude  distortion 


ARG  F(jft) 


ideal 


Q 


b)  Phase  angle  distortion 


Fig.  3.2.   FORWARD  EULER  FREQUENCY 
DISTORTION. 
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or  equivalently 


/24 
B    «   -4?~  (3-26) 


/2~4  ~    3 

ft    <<   ^4^  uis    ~    t  a)      .  (3-27) 

2tt  4      s 


The  meaning  of  Eq.  (3-27)  can  be  made  more  explicit.   Suppose 
there  exists  a  continuous  frequency  spectrum  whose  highest 
frequency  component  is  to,  .   In  order  for  the  Euler  integration 
scheme  to  transform  this  spectrum  into  the  discrete  frequency 
domain  without  distortion,  the  sampling  frequency  (to  )  must 
be  much  larger  than  4/3  of  the  corresponding  highest  digital 
frequency  component  (ft,  )  .   If  this  condition  is  satisfied  then 
the  transformation  will  be  linear.   This  implies  that  if  a 
series  of  integrators  are  used  in  the  continuous  representation, 
then  each  transformation  of  these  integrations  to  the  digital 
domain  must  satisfy  the  criterion,  namely, 

ws  >>  |  ftR  (3-28) 

Since  T  =  2ir/to  ,  Eq.  (3-28)  puts  a  constraint  on  the  value  of 
the  integration  time  step  if  a  linear  transformation  is  desired. 
However,  there  is  a  much  more  fundamental  restriction  which 
must  be  placed  on  T  as  the  following  discussion  will  show. 

The  use  of  the  Euler  transformation  should  be  limited  to 
those  continuous  filters  whose  frequency  spectrum  is  essentially 
band-limited.   This  is  a  very  important  result.   Observe  Fig. 
3.3  which  portrays  the  effect  of  the  Euler  integration  scheme 
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Fig.  3.3.   EULER  TRANSFORMATION  OF  CONTINUOUS 
FREQUENCY  SPECTRUM  TO  DISCRETE 
SPECTRUM. 
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on  transforming  a  continuous  filter  spectrum  into  the  digital 
filter  spectrum.   It  can  be  seen  that  unless  the  continuous 
spectrum  is  band-limited  the  transformation  will  cause  the 
high  frequency  ends  of  the  spectrum  to  be  "chopped-of f "  when 
being  transformed  to  the  digital  representation.   Since  the 
transforming  sine  wave  has  finite  amplitude  the  continuous 
signal  must  have  finite  bandwidth.   In  addition  the  amplitude 
of  the  sine  wave  is  directly  related  to  this  finite  bandwidth 
in  a  manner  analogous  to  the  Sampling  Theorem  [13] . 

This  theorem  states  that  for  a  signal  of  finite  bandwidth, 
W,  radians  per  second,  a  sampling  rate  equal  to  at  least  twice 
the  maximum  frequency  (W)  of  the  signal  is  necessary  to  recover 
the  signal  without  distortion  from  an  ideal  low-pass  filter. 
Fig.  3.4  represents  this  process.   Fig.  3.4a  shows  a  band- 
limited  signal  of  W,  radians  per  second.   If  the  signal  is  • 
sampled  at  a  rate  higher  than  that  specified  by  the  theorem, 
Fig.  3,4b  is  the  result.   However,  if  the  signal  is  sampled 
at  a  rate  less  than  the  Nyquist  rate,  "aliasing"  or  "folding" 
occurs  and  this  is  illustrated  in  Fig.  3.4c. 

To  relate  this  to  the  Euler  method,  consider  again  Fig.  3.3. 
If  the  continuous  filter  has  as  its  highest  frequency  component, 
go,  ,  then  to  ensure  that  no  part  of  this  spectrum  will  be  lost 
during  the  transformation,  the  value  of  T  must  be  such  that 

a)   >  2  a).  (3-29) 

s  —    h 

As  long  as  this  condition  is  satisfied  the  amplitude  of  the 
transforming  sine  wave  will  be  sufficiently  large  to  guarantee 
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b)  Spectrum  of  sampled  band 
limited  signal  2frf   >  2W, 
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c)  Spectrum  of  sampled  hand 
limited  signal,  2?rf   <  2W 


Fig.  3.4.   EFFECT  OF  SAMPLING  FREQUENCY  ON 
THE  FREQUENCY  SPECTRUM  OF 
SAMPLED  SIGNAL. 
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that  all  of  the  continuous  spectrum  will  be  mapped  into  the 
digital  domain.   This  then  is  the  fundamental  constraint  that 
is  placed  on  the  value  of  T.   This  must  always  be  satisfied 
when  using  the  Euler  method.   However,  it  is  necessary  to  set 
the  value  of  T  lower  than  that  implied  by  Eq .  (3-29)  in  order 
to  keep  the  distortion  reduced. 

Since  the  Euler  method  must  be  restricted  to  those  contin- 
uous filters  which  have  a  band-limited  spectrum,  it  is  obvious 
that  this  technique  is  limited. 

B.   BACKWARD  EULER  INTEGRATION  METHOD 

The  backward  Euler,  as  was  pointed  out  in  the  previous 
section,  differs  from  the  forward  Euler  algorithm  in  the  way 
the  input  is  treated.   The  backward  method  uses  the  present 
value  of  the  input  to  compute  the  present  value  of  the  output. 
Recall  Eq.  (3-3) 


Yk  =  Yk-1  +  T  f (Tk)  (3-30) 


and  it  can  be  seen  that  the  input  and  output  are  considered  at 
the  same  time  point.   Intuitively  this  would  appear  to  be  more 
accurate  and  would  not  involve  any  prediction  on  the  part  of 
the  integration  scheme,  and  in  general  this  is  true. 

An  alternative  way  to  understand  the  difference  between  the 
two  Euler  methods  is  to  look  at  the  differential  equation 

y  =  fit]  (3-31) 
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The  Euler  method  integrates  this  type  of  differential  equation 
by  approximating  the  derivative  as  shown 

•    Yk  ""  Yk-1 

y  =     T  K  ±  (3-32) 

In  the  forward  method  the  derivative  is  approximated  at  time 
T,  ,  and  the  formula  of  Eq.  (3-4)  results.   But  the  backward 
Euler  approximates  the  derivative  at  time  T,  and  Eq.  (3-3)  is 
applicable.   When  selecting  an  algorithm  it  must  be  decided 
whether  or  not  the  present  value  of  the  input  will  be  available 
to  compute  the  present  value  of  the  output.   It  would  be  pre- 
ferable to  use  the  present  value  of  the  input,  but  there  can 
be  circumstances  which  make  this  impossible. 

Similar  to  the  procedure  followed  in  the  previous  section 
the  z-transform  can  be  applied  to  Eq.  (3-30)  and  the  result, 
is 

Y(z)  =  z"1  Y(z)  +  T  F(z)  (3-33) 

From  this  equation  the  integration  transfer  function  can  be 
formed. 

|  (z)  =  — ~I  (3-34) 

1-z 

Comparing  this  to  Eq.  (3-11)  yields 

1-z""1 
s  -  ±-| (3-35) 

for  the  backward  Euler  integration  algorithm.   To  find  the 
sinusoidal  steady-state  response  of  Eq .  (3-35)  the 
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substitutions  mentioned  previously  are  applied.   This  results 
in  the  following 

F(jfl)  =  i=^= (3-36) 


F(jJ2)  =  i  [sin  ftT  -  j (1-cos  OT)]  (3-37) 

Comparing  this  expression  to  Eq.  (3-15),  it  is  seen  that  only 
the  sign  of  the  real  part  has  changed.   This  does  not  affect 
the  magnitude  of  F(jfi)  which  is  the  same  as  Eq.  (3-16).   Thus 
all  of  the  analysis  that  was  presented  concerning  the  frequency 
distortion  for  the  forward  Euler  is  true  for  the  backward 
Euler  and  the  same  constraints  on  the  value  of  the  integration 
time  step  are  applicable. 

The  phase  angle  characteristics  of  the  backward  case, 
however,  are  different  from  the  forward  Euler  and  are  given 
by 

ARG  F(jft)  =  j   -  ~  (3-38) 

This  is  plotted  in  Fig.  3.5  along  with  the  magnitude  distor- 
tion.  It  can  be  seen  that  the  phase  angle  of  the  backward 
Euler  is  approximately  the  same  as  the  ideal  case  when 

TT      fiT 

y  >>  -y-    .   Thus  for  low-pass  filters  with  low  cutoff  frequencies 
the  phase  angle  of  the  backward  Euler  is  virtually  the  same  as 
the  ideal  integration. 

One  property  that  would  be  desirable  for  any  integration 
scheme  to  possess  would  be  the  mapping  of  the  imaginary  axis 
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b)  Phase  angle  distortion. 


Fig.  3.5.   BACKWARD  EULER  FREQUENCY 
DISTORTION. 
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of  the  s-plane  on  to  the  unit  circle  of  the  z-plane.   The  unit 
circle  in  the  z-plane  establishes  the  region  of  stability  for 
the  Z-domain;  namely,  the  interior  of  this  unit  circle.   If  an 
algorithm  mapped  the  imaginary  axis  of  the  s-plane  onto  the 
unit  circle  of  the  Z-plane,  then  it  would  be  possible  to 
ensure  that  as  long  as  the  continuous  filter  was  stable  then 
the  digital  representation  would  also  be  stable. 

To  examine  this  characteristic  closer  consider  the  back- 
ward Euler  integration  transfer  function  Eq.  (3-35) .   Multiplying 
numerator  and  denominator  by  z  yields 

s  =  i  5§i  (3-39) 

Solving  for  z 

z  =  i=if  '  (3-40) 

Evaluating  this  expression  with   s  =  joo 


(3-41) 


1-jwT 

it  is  obvious  that  the  magnitude  of  z  does  not  equal  unity. 
In  fact 


21  ■  -   2,2,1/2  (3"42) 

(1+w  T  )  / 


ARG  z  =  tan  -1  ooT  (3-43) 


The  above  mapping  pair  transforms  the  imaginary  axis  of  the 
s-plane  to  a  pair  of  circles  that  are  completely  enclosed  by 
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b)  Backward  Euler  mapping. 

Fig.  3.6.   EULER  MAPPING  OF  Z-PLANE 

IMAGINARY  AXIS  ONTO  z -PLANE. 
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the  unit  circle.   This  is  plotted  in  Fig.  3.6b.   This  ensures 
that  the  backward  Euler  will  map  continuous  system  poles 
located  in  the  left  half  of  the  s-plane  into  the  stable  region 
of  the  z-plane.   The  forward  Euler,  however,  does  not  perform 
in  this  manner  as  shown  in  Fig.  3.6a.   The  appropriate  equations 
for  the  forward  scheme  are 


z|  =  (1  +  co2T2)1/2  (3-44) 


ARG  z  =  tan  -1  coT  (3-45) 


It  is  seen  that  the  imaginary  axis  is  mapped  onto  a  curve  which 
goes  to  infinity  along  asymptotes  of  +  90  degrees.   Thus  with 
the  forward  Euler  algorithm  there  is  no  assurance  that  a  pole 
in  the  left  half  of  the  s-plane  will  be  mapped  into  the  stable 
region. 

The  foregoing  analysis  would  seem  to  indicate  that  the 
backward  Euler  integration  would  be  superior  to  the  forward 
algorithm  for  use  in  the  algebraic  substitution  method. 
Although  they  possess  the  same  frequency  magnitude  characteris- 
tics the  phase  angle  characteristics  of  the  backward  Euler  are 
superior  to  that  of  the  forward.   Certainly  as  far  as  stability 
is  concerned  the  backward  scheme  is  superior.   The  next  sec- 
tion discusses  the  trapezoidal  integration  algorithm. 

C.   TRAPEZOIDAL  INTEGRATION  METHOD 

The  trapezoidal  method  of  discrete  integration  derives  its 
name  from  the  geometric  shape  of  the  approximation  that  is 
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made  to  the  function  being  integrated.   Instead  of  using  a 
piecewise  constant  approximation  (rectangles)  as  was  done  by 
the  two  Euler  methods,  the  trapezoidal  scheme  uses  a  trapezoid 
to  fit  the  function  as  shown  in  Fig.  3.7.   This  results  in  a 
more  accurate  fit  to  the  actual  function  and  reduces  the  error 
between  the  actual  and  computed  values.   Recalling  Eq.  (3-5) 


\  -  Yk-i + !  i«v  +  f^k-i>]  (3-46) 


it  can  be  seen  that  for  the  trapezoid  algorithm  the  solution 
to  the  differential  equation 

y  =  f(t)  (3.47) 

is  being  approximated  at  the  average  of  the  two  time  points 
T,  and  T,  - .   This  characteristic  ensures  a  "good  fit"  for 
functions  which  are  montonically  increasing  or  decreasing. 
Also  this  algorithm  is  quite  good  if  the  input  function  is 
reasonably  well-behaved. 

To  obtain  the  algebraic  substitution  for  the  trapezoidal 
method  it  is  necessary  to  z-transform  Eq.  (3-46)  and  doing 
this  yields 

Y(z)  =  z"1  Y(z)  +  |  [F(z)  +  z"1  F(z)]  (3.48) 

The  integration  transfer  function  can  now  be  formed  and  is 

|  Cz)  =  |±±£i  (3-49) 

1-z 
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'f  (t) 


i t * 


T, 


y  = 


k-l 


f(t)  dt 


k-l 


Yk  =  Yk-1  +  2    [f(Tk}  +  f(Tk-l)] 


Fig.  3.7.   TRAPEZOIDAL  APPROXIMATION  TO 
AREA  UNDER  CONTINUOUS  CURVE. 
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An  examination  of  Eq.  C3-49)  reveals  that  the  integration 
transfer  function  possesses  a  bilinear  nature.   In  fact  several 
authors  [1,  5]  refer  to  the  trapezoidal  transformation  as  the 
bilinear  substitution.   This  bilinearity  is  one  property  of 
this  algorithm  which  makes  it  very  useful  for  digital  filter 
design.   Proceeding  as  before,  it  is  possible  to  make  the 
correspondence  between  the  s-plane  and  the  z-plane 

2   1-z"1 
s  *  |  i-^j.  (3-50) 

1-z 

This  is  the  algebraic  substitution  for  the  trapezoidal  inte- 
gration scheme. 

By  algebraic  manipulation  of  Eq .  (3-50),  one  of  the  most 
desirable  features  of  the  bilinear  transformation  can  be 
demonstrated.   It  was  stated  earlier  that  it  would  be  convenient 
for  a  transformation  to  map  the  imaginary  axis  of  the  s-plane 
onto  the  unit  circle  of  the  z-plane.   The  trapezoidal  algorithm 
possesses  this  property.   To  see  this  consider  Eq.  (3-50)  and 
disregard  the  constant  factor  2/T. 

z-1 


"    z+1 
Solving  for  z  yields 


(3-51) 


z  =  i±5.  (3-52) 

1-s 


Let  s  =  joo 


z  .  *£j»    +        |z|  -  1  (3-53) 

1-jo)        ■  ■ 
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This  transformation  is  the  only  one  which  possesses  this 
property  of  mapping  the  imaginary  axis  onto  the  unit  circle 
without  increasing  the  order  of  the  system.   This  ensures  that 
a  pole  in  the  left  half  of  the  s-plane  in  the  continuous 
system  will  be  mapped  into  the  stable  region  of  the  z-plane. 
This  fact  makes  the  trapezoidal  integration  algorithm  a  very 
useful  tool  for  the  algebraic  substitution  method. 

To  examine  the  frequency  domain  characteristics  of  the 
trapezoidal  method,  apply  the  substitutions  mentioned  in 
Section  A.  to  Eq.  (3-50) .   The  result  is 

2    l-e"jnT 

F(^>   mT7f=W  (3"54) 

1+e    J 


F(j'fi)    =    j|     tan   ^  (3-55) 


The  first  item  of  interest  to  be  noted  is  that  the  F(jft)  has 
no  real  part  and,  consequently,  the  phase  of  the  trapezoidal 
algorithm  is  precisely  the  same  as  that  of  the  ideal  operator, 
This  characteristic  is  very  important  when  developing  the 
frequency  substitutions  as  will  be  seen  in  Chapter  IV.   To 
see  the  magnitude  distortion  compare  Eqs .  (3-19)  and  (3-55) 

j«  -  j  £  tan  f-  (3-56) 

^  =  tan  01  (3-57) 

2         2 

This  magnitude  distortion  is  clearly  shown  in  Fig.  3.8.   It 
can  be  seen  from  Fig.  3.8  that  in  the  low  frequency  case  the 
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T" 
TT 

T 


ft 


Fig.  3.8.   TRAPEZOIDAL  ALGORITHM  FREQUENCY 
DISTORTION. 
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transformation  function  can  be  considered  linear  and,  in  fact, 
this  transformation  is  linear  for  values  of  ft  such  that 

ft  <<  j   o>s  (3-58) 

The  above  condition  is  similar  to  that  given  for  the  Euler 
transformations  in  the  previous  section.   Both  conditions  result 
from  the  fact  that  the  trigonometric  terms  involved  are  linear 
about  the  origin.   However  there  is  an  important  difference 
between  the  trapezoidal  and  the  Euler  methods  which  is  also 
attributable  to  these  trigonometric  terms.   Fig.  3.9  shows 
the  transformation  from  a  continuous  spectrum  to  a  digital 
spectrum.   The  fact  that  the  trapezoidal  transformation 
contains  the  tangent  factor  means  that  the  entire  continuous 
frequency  band  can  be  transformed.   In  fact,  the  entire 
continuous  spectrum  is  transformed  into  the  interval  +  tt/T 
in  the  discrete  domain.   This  property  makes  the  bilinear 
substitution  a  more  universal  type  of  transformation  as  it 
can  be  used  for  any  spectrum. 

Unfortunately,  for  the  filter  designer,  there  still  exists 
the  problem  of  distortion  when  using  the  trapezoidal  method. 
This  distortion  can  be  overcome,  in  part,  by  pre-warping  [5]  or 
frequency  scaling  the  particular  continuous  realization  so  that 
after  the  algebraic  substitution  is  made  the  selected  critical 
frequencies  are  at  the  desired  values.   However  this  will  work 
only  for  those  frequencies  selected  and  may  yield  an  unsatis- 
factory design  at  other  frequencies.   Another  solution  and  one 
that  has  been  mentioned  previously  would  be  to  increase  the 
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Fig.  3.9.   TRAPEZOIDAL  TRANSFORMATION  OF  CONTINUOUS 
FREQUENCY  SPECTRUM  INTO  DISCRETE  DOMAIN. 
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sampling  frequency ,  thereby  ensuring  that  the  critical 
frequencies  would  be  much  less  than  the  sampling  frequency. 
This  solution  also  has  drawbacks.   Increasing  the  sampling 
frequencies  also  reduces  the  time  available  for  computation. 
If  the  designer  is  faced  with  a  fixed  computing  speed  then 
increasing  the  sampling  frequency  would  limit  the  complexity 
of  the  filter  than  can  be  implemented.   Thus  reducing  the 
sampling  interval  might  not  be  a  desirable  tool  to  use.   In 
this  case  pre-warping  might  be  the  only  feasible  solution. 

The  trapezoidal  or  bilinear  transformation  has  been  used 
more  extensively  than  any  other  substitution  technique  due  to 
its  universal  application  to  the  frequency  domain  and  the  fact 
that  it  does  not  increase  the  complexity  of  the  filter;  i.e., 
an  n    order  continuous  filter  is  transformed  into  an  n    order 
discrete  filter.   Also  the  fact  that  the  trapezoidal  algorithm 
possesses  ideal  phase  angle  characteristics  greatly  simplifies 
the  frequency  substitution  as  will  be  demonstrated  in  Chapter 
IV. 
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IV.   NUMERICAL  INTEGRATION  - 
GENERAL  ALGORITHM 

Chapter  III  presented  the  general  case  discrete  integration 
formula  and  it  is  repeated  here  for  emphasis 


Yk  =  Yk-i  +  T  x  f  (Tk-i)  +  T^-X^  f  <V  (4_1) 


The  two  quantities  in  this  formula  that  are  of  particular 
interest  are  T  (the  integration  time  step)  and  A.   The  parameter 
A  may  be  considered  a  weighting  parameter  in  that  the  selection 
of  X  determines  which  time  point  of  f(t)  will  have  the  greatest 
effect  on  the  integration  process.   These  two  quantities  (T 
and  A)  determine   the  type  of  performance  that  can  be  achieved 
using  this  integration  formula.   It  was  seen  earlier  that 
choosing  the  values  of  A  to  be  0,  ~-,  and  1  led  to  three  well- 
known  integration  algorithms.   However  there  is  no  absolute 
restriction  on  this  value  of  A  and  there  may  be  instances 
where  another  choice  would  provide  better  results. 

It  is  known,  however,  that  a  value  of  A  such  that 

0  1  x  1   |  (4"2) 

guarantees  stability  of  the  numerical  integration  process. 
This  constraint  on  the  value  of  A  means  that  as  long  as  the 
poles  of  the  continuous  system  lie  in  the  left  half  of  the 
s-plane,  the  solution  to  the  system  will  be  numerically 
stable  for  any  positive  value  of  T.   If  A  >  y  then  the  value 
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of  T  will  have  some  upper  bound  that  will  depend  on  the  exact 
value  of  A  that  is  used  and  on  the  location  of  the  roots  of  the 
system.   Appendix  B  contains  a  detailed  analysis  of  these 
stability  conditions.   It  is  sufficient  here  to  be  aware  of 
the  advantages  that  accrue  from  choosing  A  according  to 
Eq.  (4-2)  . 

Fig.  4.1  is  a  representation  of  the  general  integration 
formula  for  various  values  of  A .   It  can  be  seen  that  the 
approximation  to  the  actual  function  moves  up  or  down  as  A 
is  varied.   Thus  as  A  is  made  smaller  the  approximating 
rectangle  grows  larger  and,  conversely,  as  A  is  made  larger 
the  approximating  rectangle  becomes  smaller.   Thus  it  is  seen 
that  the  value  of  A  determines  which  time  point  will  have  the 
greatest  effect  on  the  approximation. 

It  would  be  convenient  to  have  the  general  integration 
formula  available  for  use  in  the  algebraic  substitution  method. 
This  would  greatly  facilitate  implementing  the  method  once  the 
designer  has  selected  the  parameters  he  wants  to  use.   This 
chapter  contains  the  derivation  of  the  general  case  transfor- 
mations that  could  then  be  employed  once  the  values  of  T  and 
A  have  been  selected.   In  addition  the  derivation  of  the 
general  case  frequency  substitutions  are  presented  in  section 
B. 
A.   GENERAL  ALGORITHM  TRANSFORMATIONS 

Following  the  procedures  outlined  previously  it  is 
necessary  to  z-transform  Eq .  (4-1)  to  put  it  in  the  form 
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k-1 


rTi 


y  = 


f(t)  dt 


k-1 


Yk  =  Yk-1  +  T  A  f(Tk-l}  +  T(1-x)  f(Tk> 


Fig.  4.1.   GENERAL  INTEGRATION  ALGORITHM 

APPROXIMATION  FOR  VARIOUS  VALUES 
OF  A. 
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required  for  the  algebraic  substitution.   Doing  this  yields 
Y(Z)  =  Z"1  Y(Z)  +  T  A  z"*1  F(Z)  +  T(l-A)  F(Z)     (4-3) 


Solving  to  form  the  integration  transfer  function  produces 

Y(Z)  [1-Z"1]  =  TUCZ"1-!)  +  1]  F(Z)  (4-4) 


|  (z)  =  TKZ-1-!)^,!] 

F        (l-z-1) 


Eq.  (4-5)  is  the  desired  integration  transfer  function  and 
from  this  the  following  correspondence  can  be  made. 

1    (l-z-1)  ,.   c, 

s  +   =•  ^ '- (4-6) 

[(Z   -1)A+1] 

This  then  is  the  required  algebraic  substitution  for  s  for  the 
general  integration  algorithm.   By  selecting  a  value  for  A 
one  can  immediately  obtain  the  desired  substitution  for  s. 
It  can  be  seen  that  when  A  =  0,  ~-,  1,  Eq.  (4-6)  becomes  in  turn 
Eq.  (3-35),  Eq.  (3-50),  and  Eq .  (3-12). 

The  next  transformation  required  is  that  for  mapping  from 
the  continuous  frequency  domain  to  the  discrete  frequency 
domain.   To  accomplish  this  consider  Eq.  (4-6)  and  let 

Z_1  =  exp  (-jftT) 
Applying  this  to  Eq.  (4-6) 

P(jfi)  =  k  -ircr (4""8) 

1  [1+A(e  ^-l)] 
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This  form  can  be  expressed  in  terms  of  trigonometric  quantities 
as  follows 

2 .        sin  2~ 

F(^}  "  T2   ^T   .  ,„  ~         .   ttT  <4~9> 

cos  y~  ~J    (2X-1)  sin  y~ 

This  form  permits  the  necessary  F(jftT)  expression  to  be  quickly 
developed  for  substitution  into  the  frequency  domain  transfor- 
mations which  are  derived  in  the  next  section. 

To  analyze  these  generalized  results  it  is  helpful  to 
rearrange  Eq .  (4-9)  into  a  form  similar  to  the  ones  used  in 
the  preceding  chapter.   Eq.  (4-9)  can  be  written 


ftT      ftT    ,.,  .,   .  2  J2T 
?  j  sin  -y—   cos  2~-  "  (2X-1)  sin   =— 

^  =  T 2  AT  .,  „,  ,.2   ■  2  BS (4"10) 

cos   2~  +  (2A.-1)   sin  -~— 

2  f^T 
Dividing  through  by  sin   ~ —  yields 


QT 
2  -(2X-1)  +  j  cot  ji- 

jw  -  ip  2^ 2"^ (4-11) 

1   (2A-1-)2  +  cot^  ~ 

Consider  the  magnitude  of  both  sides  of  Eq.  (4-11) . 

a)  =  =■  2—^ 2  1/2  (4-12) 

T  [cot:  f-  +    (2X-1)^]1/Z 


Eq.  (4-12)  permits  an  examination  of  the  frequency  distortion 
that  will  be  present  when  this  general  integration  algorithm 
is  used. 

Fig.  4.2  is  a  representation  of  Eq.  (4-12)  for  various 
values  of  X.   It  can  be  seen  that  when  X  =  0,  j,    and  1  the 

53 


2.0 


1.0 


1 

1 

1     1 

" 

1 

hi 

_LLi_i 

_. .  i     L. 

. 

! 

UJ 

._[_'    .    . 

M 

.          ._ 

1 

! 

'  1 

r 

1 

: 

. 

I 

/ 

• 

I  <  1 

i     '' 

i   1      1 

_J_L 

1 

i 

/ 

i 

" 

i 

■ 

/ 

i  ! 

■"[ 

/ 

1   ! 

1  1 

III 

! 

'   ' 

II 

j   i 

_J_L.. 

/ 

i 

! 

i 

1    1 

1 

1 

i—  i  . 

/ 

1    !    I 

; 

A=* 

/ 

! 

I 

!_' 

1 

2     V 

, . i  i 

/ 

Ml 

<M 

1 

Ml 

| 

j 

. 

II 

. 

!    1 

>v 

: 
; 

^r 

^         1 

x^ 

1 

i 

L      . 

j  M 

j 

i  i 

;  i  - 

i   !   | 

i    1    1 

1 

!     1    ! 

I 

s^ 

1      ! 

I    ! 

mTT 

1    > 

i 

f 

s\ 

I 

M~     - 

1    1 

1        1 

i  / 

V 

1 

, 

/ 

1    ! 

! 

. 

t= 

1 

.        1 

/ 

_Y 

1 

i      ! 

1 

1 

1  ' 

A 

>^ 

/ 

1    1    1 

1 

^s^ 

/ 

■ 

s^ 

' 

....:.:    : 

1 

i 

^S^ 

. 

I  j 

. 

^* 

<w/ 

. 

1 

!    ■ 

II 

. 

jj    1    I 

i  ._ 

1     1 

' 

^ffr 

/  :  y 

,     -      1 

~r 

1 

/       / 

1 

. 

1 

,    ,        I 

1 

1 

S" 

/ 

I 

.... 

l  l 

i 

/ 

/ 1        -- 

—     1 

Ml 

rill 

!  < 

1 

/  / 

/ 

i 

1 

Ml 

1        1    i 

\  1 

!       !    ! 

/  a^     1 

?  ***^ 

J 

i 

I 

> 

/ jr 

^i 

i  i  j  i 

I'll 

i    i 

/ 

, 

XI 

i  1  i 

i  !  • 

1 

1  1  i 

i  i  | 

\ 

•    : 

i  1  1  -1 

1    1       1 

i 

f 

'    1 

! 

_•    |    lJ_ 

: 

\=0    ' 

1  1    . 

lii 

£>^                                          / 

i  i  i  i  j 

^T 

1 

/ 

\  i 

:    j   ' 

i  | 

^     i 

' 

A 

:  i 

III 

i 

„_L 

1    1   I 

; 

\- 

.     .  8 

!     • 

L 

i  ■  1 

>  rr 

.  ■_. 

;   U- 

. 

! 

1  ,  1  i 

.    1  _ 



: !  1 

-jr 

. 

Ml 

TT 

2T 


2L 

T 


n 


Fig.  4.2.   PLOT  OF  FREQUENCY  TRANSFORMATIONS  FOR 
VARIOUS  VALUES  OF  A . 
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plots  are  the  same  as  in  the  previous  chapter.   However,  for 
other  values  of  X   there  are  no  clear-cut  analytical  expres- 
sions for  the  curves .   The  amplitude  of  the  curves  increase 
as  X  approaches  j   so  that  the  trapezoidal  rule  (X  =  j)  is  the 
limiting  upper  bound.   Conversely,  the  amplitude  of  the  curves 
decrease  as  X  approaches  either  0  or  1  so  that  the  Euler 
methods  are  the  limiting  lower  bound. 

In  the  previous  chapter  it  was  discovered  that  the  Euler 
transformations  should  be  used  only  for  band-limited  spectra. 
Also  the  value  of  T  that  was  used  would  have  to  be  constrained 
in  such  a  manner  as  to  ensure  that  all  of  the  spectrum  would 
be  transformed  without  loss.   Here  in  the  general  case  it  can 
be  seen  that  the  same  situation  exists.   All  of  the  transfor- 
mations except  for  the  trapezoidal  one  have  finite  amplitudes 
and  thus  should  only  be  used  for  finite  spectra.   However, 
even  though  the  amplitudes  are  finite  they  can  be  made  arbi- 
trarily large  by  the  proper  choice  of  X.   This  means  that  the 
constraint  on  the  value  of  T  may  be  relaxed.   Unfortunately, 
linearity  considerations  must  also  be  taken  into  account  and. 
these  will  cause  the  value  of  T  to  be  kept  small  to  avoid 
excessive  distortion. 

Another  factor  that  must  be  considered  at  this  point  is 
just  how  complex  a  transformation  one  would  want  to  use.   It 
is  obvious  from  Eq.  C4-6)  that  a  value  of  X  other  than  0,  y, 
and  1  will  make  the  algebra  associated  with  that  particular 
substitution  more  cumbersome.   In  addition  any  value  of  X 
other  than  ~-  will  cause  some  phase  distortion  and  this  may  be 
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an  important  consideration.   In  fact,  section  B  will  show  that 
when  the  phase  is  not  ideal,  the  frequency  transformations 
become  much  more  involved. 

B.   GENERAL  CASE  FREQUENCY  TRANSFORMATION 

It  was  pointed  out  in  Chapter  II  that  it  would  be  convenient 
to  have  a  substitution 

a)  =  F(jfi)  (4-13) 

that  could  be  used  to  directly  obtain  the  magnitude  squared 
representation  for  the  discrete  transfer  function.   For  each 
integration  algorithm  that  has  been  considered  a  F(jfi)  has 
been  derived.   Now  if  this  expression  was  the  proper  substitu- 
tion for  co  then  this  problem  would  not  only  be  solved,  but 
greatly  simplified  as  these  expressions  were  not  very  complex. 
However,  in  general,  these  expressions  do  not  fulfill  the 
task.   To  clearly  demonstrate  this,  a  simple  example  will  be 
presented. 

Consider  the  transfer  function  depicted  in  Fig.  2.2 


Jffit  (.,  -  fi^|  (4-14, 

in 


If  the  backward  Euler  transformation  is  applied  to  this 
continuous  expression,  a  transfer  function  in  the  z-domain 
is  obtained. 

fout  (z)  m   1  ±   a  T  +  z"^  (4-15) 

Ein        1  +  b  T  +  z 
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Now  let  z  =  eJ    and  obtain  the  magnitude  squared  function. 


out 


in 


(jfl) 


2       2 

CI  +  a  T  -  cos  AT)   +  sin   QT 

2        2 
CI  +  b  T  -  cos  fiT)   +  sin   fiT 


(4-16) 


This  is  the  correct  expression  as  it  was  derived  directly  from 
the  discrete  transfer  function.   However  this  procedure  would 
usually  be  very  time-consuming  and  a  better  solution  would  be 
to  directly  substitute  an  expression  into  the  magnitude  squared 
sinusoidal  steady-state  transfer  function. 


out 

e . 
in 


Cjoo) 


2     2 

oj  +  a 

2  _L  V2 

w  +  b 


(4-17) 


However  when  this  is  done  using  the  F(jfi)  derived  earlier  the 
function  obtained  is 


E 


out 


in 


(jfl) 


2  (1  -  cos  AT)  +  a2T2 
2  (1  -  cos  QT)    +   b2T2 


(4-18) 


It  is  evident  after  only  a  little  algebra  that  Eq .  (4-16)  and 
Eq .  (4-18)  are  not  equivalent,  and  consequently  the  procedure 
employed  is  not  valid  in  this  case. 

The  negative  result  of  the  above  example  demonstrates  the 
need  for  a  general  transformation  which  will  be  valid  for 
every  case.   Two  separate  cases  were  considered  in  developing 
these  general  transformations.   The  first  case  was  the  linear 

term  (s  +  a)  as  was  used  in  the  preceding  example.   The  second 

2 

was  the  quadratic  factor  (s   +  as  +  b) .   Since  transfer  function 
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factors  should  be  put  into  no  more  than  second  order  terms 

for  ease  in  implementing  the  filter  it  was  considered  sufficient 

to  derive  general  transformations  for  these  two  cases  only. 
Consider  the  linear  term  situation. 

(s  +  a)  =  [fCz)  +  a]l  =  f(z)  (4-19) 

Now  examine  the  magnitude  squared  function. 

w2  +  a2  =  |F(ejfiT)  +  a|2  f(z)  =  F(eJOT}  (4-20) 

=  (a  +  Re  F(ejfiT))2  +  lm2  F(ejOT)      (4-21) 
w2  +  a2  =  a2  +  2a  Re  F(ejfiT)  +  |F(e^T)|2       (4-22) 


CO 


2  =  2a  Re  F(ejfiT)  +  |F(ejfiT)|2  (4-23) 


This  then  is  the  final  result  for  the  linear  case.   It  can  be 
seen  that  the  term  (2a  Re  F(eJ   ))  is  the  one  which  prevents 
the  simpler  transformation 

2 

2     i     -iQT  I 
oo   =   Fte-5"1)  (4-24) 


from  being  valid  in  the  general  case.  For  those  transforma- 
tions which  do  not  have  a  real  part,  e.g.,  trapezoidal  rule, 
the  general  transformation  of  Eq.  (4-23)  reduces  to  the  much 
simpler  form  of  Eq.  (4-24) .  It  can  also  be  seen  that  when  a 
transformation  does,  have  a  real  part,  this  is  equivalent  to 
saying  that  the  phase  angle  characteristics  differ  from  those 
of  the  ideal  integrator.   This  is  the  case  with  both  Euler 

methods . 
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The  case  of  the  quadratic  factor  requires  considerable 
mathematical  manipulation  and  yields  a  transformation  which  is 
quite  a  bit  more  complex.   Assuming  a  quadratic  factor  of  the 
form 


F(s)  =  s2  +  a  s  +  b  (4-25) 


the  magnitude  squared  function  becomes 


F(tD2)  =  u)4  +  oo2  (a2  -  2b)  +  b2  (4-26) 


By  a  procedure  similar  to  that  followed  for  the  linear  case 
the  following  transformation  was  obtained 


4     2,  2  0,x    ,2    i_2,  jfiT.  i2     2  ._,  jOT.,2,,2 
u   +  ca  (a  -2b)  +  b   =   F  (eJ   )    +  a    F(eJ   )   +b 


+  2a  Re  F3(ejfiT)  +  2b  R  F2  (ejfiT) 


+  2a  b  Re  F(e^T)  (4-27) 


Thus  the  entire  factor  in  the  continuous  domain  becomes  that 
shown  in  Eq.  (4-27)  in  the  discrete  domain.   Even  in  this  case, 
however,  for  those  algorithms  which  possess  ideal  phase 
characteristics  the  general  transformation  reduces  to  a  much 
simpler  and  more  logical  form.   Unfortunately,  in  general,  a 
simple  substitution  will  not  suffice. 

These  transformations  provide  a  direct  avenue  to  forming 
the  frequency  domain  representations  of  a  digital  filter. 
While  they  do  not  appear  simple  in  form,  in  many  cases  they  will 
reduce  to  a  very  convenient  expression  which  will  provide  the 
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designer  with,  the  desired  information  regarding  the  frequency 
response  of  a  particular  filter.   They  will  also  provide  a 
good  deal  of  insight  into  the  frequency  distortion  that  will 
be  present.   This  is  demonstrated  in  the  next  section. 

C.   FREQUENCY  DISTORTION  CONSIDERATIONS 

It  has  been  shown  that  every  algorithm  that  can  be  derived 
from  the  general  integration  formula  possesses  some  degree  of 
frequency  distortion.   In  Chapter  III  this  distortion  was 
explicitly  analyzed  for  three  specific  algorithms.   In  addi- 
tion guidelines  were  presented  for  minimizing  this  distortion. 
These  rules  indicated  that  by  decreasing  the  sampling  interval 
(T)  appropriately,  the  amount  of  distortion  can  be  reduced  to 
a  negligible  value  over  a  specified  frequency  range.   In 
general  this  is  the  case  with  any  algorithm  that  is  used.   As 
long  as  the  sampling  frequency  is  much  greater  than  the  cut- 
off frequency  the  transformation  involved  will  be  essentially 
linear.   However,  the  general  integration  formula  also  involves 
the  parameter  X    and  certainly  it  would  be  helpful  to  know  which 
value  of  X  is  optimum  for  reducing  distortion. 

To  examine  this  consideration  further  recall  Eq .  (4-12) 


u  -  £  0-7^ *r~- r^s-  (4-28) 


T   r   .2  ftT    ,  ..  ,.2,1/2 
[cot  ~—   +  (3X-I)  ]  ' 


This  equation  can  also  be  written  as 

ftT 
0        tan  ~— 

u  =  I  2 2  QT  1/2  (4"29) 

T  [1  +  (2A-1)2  tan2  Si]  ' 
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At  this  point  it  is  necessary  to  apply  series  expansion 
approximations .   These  approximations  are  valid  as  the  function 
here  is  being  examined  in  the  vicinity  of  the  origin  and  the 
argument  (— )  is  very  small.   Hence  higher  order  terms  can  be 
discounted.   This  leads  to 

(tiT         (ftT)3 

to  =  I  i 2 T (4-30) 

[1  +  (2X-1)2(^+  (§S?  )2]V2 


Then 


,,  *  2  r^T    (fiT)3iri    1  ,-,  ...2,fiT    (QT)3  2,  ,.  _., 
T  ^2~~  +   2~~  ^       ~    2     ^2A"1^  ^2~~    6~   *  ^  (4-31) 


since   ytt   ~  1  ~  o   provided   a  <<  1 

(l+a)1/z        ^ 


ui  ~  |  [01  +  (-|  (2X-1)2  +  |)  (^)3]  (4-32) 


At  this  point  it  is  possible  to  form  the  ratio  of  the  nonlinear 
term  and  the  linear  term  to  obtain  a  representation  of  the 
fractional  deviation  from  linearity. 

Fractional  Deviation  =  A  ■  [j   -  |  (2X-1)2]  (~)  2     (4-33) 

The  item  of  interest  that  can  be  seen  from  Eq .  (4-33)  is  that 

the  value  of  A  will  be  small  since  it  involves  a  second  order 

11     2 
term.   Also  the  factor  [~-  -  yC2A-l  )]  can  be  analyzed  for 

different  values  of  X  and  the  result  is  very  interesting.   The 

value  of  X  that  gives  the  largest  first  order  distortion 
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is  j   and  there  are  two  values  that  result  in  no  first  order 
distortion.   These  are  0.09  and  0.91.   Obviously  one  would  not 
be  interested  in  the  larger  value  since  it  would  not  insure  a 
stable  integration.   However,  the  smaller  value  would  insure 
this  stability  while  not  contributing  any  first  order 
distortion. 

While  this  analysis  is  only  valid  in  the  vicinity  of  the 
origin  it  must  be  remembered  that  the  designer  will  want  to 
adjust  the  sampling  interval  so  that  the  approximation  v/ill  be 
in  this  linear  zone.   The  important  point  of  this  analysis  is 
that  the  designer  will  be  faced  with  decisions  about  which 
algorithm  to  use  and  it  may  be  necessary  to  make  allowances 
for  a  certain  amount  of  distortion  in  order  to  have  a  less 
complex  transformation.   However,  the  analysis  that  has  been 
presented  will  permit  the  designer  to  make  intelligent  choices 
concerning  the  necessary  parameters  in  employing  the  algorithms 
At  this  point  it  would  be  instructive  to  look  at  an  example 
of  a  filter  design  problem  to  gain  insight  into  how  the  theory 
applies  in  practice. 

Example .   The  continuous  transfer  function  mentioned  in  the 
previous  section  can  be  used  to  understand  the  method  of 
attacking  the  approximation  problem.   The  magnitude  squared 
function  of  that  filter  was 

2     2 

H(w2)  =  a)2  +  a2   (a  <  b)  (4-34) 

(o  +  b 
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Also  the  discrete  counterpart  can  be  formed  using  the  transfor- 
mation derived  in  section  B. 


2       2 
r,/-.0\|2    (1  +  a  T  -  cos  ftT)   +  sin   ftT        ,.  __, 

GC3^)|   =  2 5 (4-35) 

(1  +  b  T  -  cos  SIT)       +  sin   ftT 


This  can  be  written  as  follows. 


4  M  mN   .  2  fiT     2 

-y  (1  +  a  T)  sin  j-  +  a 

G(JR)  '2  =  4  n  ^  „  „.   ■  2  W  I  h2            <4"36> 

— k-  (1  +  b  T)  sin  y-  +  b 

T^  ^ 


Fig.  4.3a  shows  a  plot  of  both  Eq .  (4.34)  and  Eq .  (4-36)  where 
Eq.  (4-36)  is  merely  shown  in  general  form  for  any  arbitrary 
value  of  T. 

Since  this  particular  continuous  spectrum  is  not  band- 
limited,  the  Backward  Euler  transformation  was  not  the  optimal 
one  to  use.   However,  its  use  illustrates  the  ideas  that  are 
involved  and  also  show  that  part  of  the  spectrum  is  lost  as 
a  result  of  the  transformation.   Fig.  4.3a  clearly  shows  the 
periodicity  of  the  filter  frequency  response.   In  addition  it 
is  interesting  to  note  that  if  T  is  allowed  to  approach  zero, 
the  digital  filter  response  approaches  that  of  the  continuous 
filter.   This  is  of  course  consistent  with  the  theory  of 
sampling.   This  also  confirms  the  assertion  made  earlier  that 
for  very  small  values  of  the  sampling  interval  the  transforma- 
tion is  essentially  linear. 

Thus  it  is  possible  to  reproduce  the  continuous  filter 
frequency  response  as  close  as  necessary  if  the  value  of  T 


63 


H(jco)  |2/|G(jr2)  I2 


w/n 


a)  Comparison  of  frequency  response  for 
digital  and  continuous  filter  without 
adjusting  sampling  internal. 


H(jco)  |2/|G(jfi)  |2 


b)  Continuous  and  discrete  spectra  with 
sampling  interval  reduced  to  achieve 
linearity . 


Fig.  4.3.   'FREQUENCY  RESPONSE  COMPARISONS  FOR 
FILTER  DESIGN  PROBLEM. 
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is  adjusted  properly.   For  this  example  if  one  wished  to  have 
the  digital  filter  response  approximate  that  of  the  continuous 
filter  for  values  of  co  up  to  ten  then  a  value  of  T  such  that 
T  =  0.01  will  do  the  job.   This  is  shown  in  a  magnified  view 
in  Fig.  4.3b.   Similarly  any  such  specification  can  be 
satisfied  if  the  designer  is  willing  to  make  the  value  of  T 
small  enough. 

The  transformations  covered  in  this  chapter  offer  the 
designer  a  variety  of  approaches  to  the  design  problem.   While 
one  transformation  may  achieve  greater  linearity  than  another, 
the  second  may  be  easier  to  implement.   Therefore,  the  filter 
designer  is  faced  with  the  recurrent  engineering  problem. 
What  selection  of  the  available  choices  can  be  put  together 
to  give  the  best  possible  design?   The  analyses  presented  will 
make  the  decisions  less  difficult. 
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V.   SUMMARY  AND  SUGGESTIONS 
FOR  FURTHER  RESEARCH 

The  important  points  presented  in  this  thesis  are: 

1)  The  algebraic  substitution  method  is  a  convenient  and 
relatively  simple  procedure  to  use  in  designing  a  digital  filter 
from  continuous  filter  characteristics.   While  the  discrete 
filter  is  only  an  approximation  to  the  continuous  one,  it  is 
possible  to  make  this  discrete  realization  approach  the 
desired  realization  as  close  as  necessary. 

2)  The  problem  of  frequency  distortion  is  present  in  any  of 
the  transformations  that  would  be  used  in  the  algebraic 
substitution  technique.   However,  the  effect  of  this  distortion 
can  be  offset  by  intelligent  application  of  the  guidelines 
contained  in  this  thesis.   These  guidelines  include  the  rules 
for  achieving  the  desired  amount  of  linearity  in  the  transforma- 
tions, the  s+ability  constraints,  and  the  general  case  frequency  transformations. 

The  foregoing  studies  indicate  that  an  interesting  area 
for  further  exploration  and  analysis  would  be  a  computer 
study  of  these  various  substitution  techniques  employing  the 
theory  presented  herein.   This  would  involve  the  use  of 
different  algorithms  with  different  classical  filters. 
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APPENDIX  A 
STABILITY  CONSIDERATIONS 

An  important  aspect  of  any  numerical  integration  method 
is  the  stability  of  the  method  [14].   In  this  context  stability 
refers  to  whether  or  not  the  errors  generated  at  each  step  of 
the  numerical  integration  process  tend  to  decay  (stable)  or 
grow  (unstable)  as  the  solution  progresses  in  time.   The 
question  of  stability  becomes  important  in  the  selection  of  how 
small  an  integration  time  step  must  be  used.   It  would  be  very 
convenient  for  an  integration  algorithm  to  be  stable  for  all 
positive  values  of  T  as  long  as  the  eigenvalues  of  the  system 
being  processed  are  in  the  left  half  of  the  s-plane. 

To  develop  a  general  stability  criterion  consider  the 
general  algebraic  substitution. 

=  1     1-z"1  (A.l) 

T   (z"1-l)X+l 

Solving  for  z. 

STX  +  1  (A. 2) 

z    STX-ST+1 

Now  substitute  a  root  at  s  =  a+jw. 

XT  a  +  XTjco  +  1  (A.  3) 

z  ~  Ta+jwT(A-l)  +  1 
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For  stability  |z|<  !•  or 


?      J     0     9 
U+^TA)       T  <  1  (A.  4, 


(1+TaX-Tcr)2  +   (ooAT-ooT)2 


After  much  algebra  the  result  is: 

T(1-2X)  >   |  a  (A. 5) 

For  a  root  in  the  left  half  of  the  s-plane  the  value  of  a  will 
be  negative.   Therefore  the  quantity  on  the  right  becomes  a 
negative  number.   The  quantity  on  the  left  will  be  greater 
than  or  equal  to  zero  as  long  as 

0  <  X    <_   j  (A. 6) 

Therefore  as  long  as  Eq.(A.6)  is  satisfied  then  any  positive 
value  of  T  will  ensure  the  inequality  in  Eq.(A.5)  is  satisfied 
and  the  integration  will  be  stable.   This  does  not  imply  that 
the  solution  will  be  unstable  for  X    >   -^    ,    but  rather  that  a 
constraint  on  the  size  of  T  must  be  imposed  to  ensure  stability 
This  is  precisely  the  case  with  the  forward  Euler.   However, 
the  Trapezoidal  and  backward  Euler  are  completely  stable  in 
this  sense. 
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